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ABSTRACT 


"  The  DFT-based  maximum  likelihood  method  (MLM)  Cor 
frequency  estimation  oC  complex  sinusoids  as  proposed  in  [1,2] 
is  extended  to  treat  the  case  of  real  data.  Improvements  on 
estimation  precision  and  computation  efficiency  are  obtained  by 
imposing  an  equal-frequency  constraint  on  the  pair  of  complex 
sinusoids  which  corresponds  to  a  real  sinusoid  and  by  utilizing 
a  single-step  interpolation  in  conjunction  with  a  coarse  finite 
search  over  the  DPT  data  for  the  maximum  of  the  likelihood 
functon.  Simulation  results  demonstrate  these  improvements. 
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I.  INTRODUCTION 

The  maximum-likelihood  (ML)  method  for  frequency 
estimation  of  sinusoids  in  noise  was  first  proposed  in  [1]  and 
extended  in  [2,3].  It  has  been  demonstrated  that  this  method 
is  often  an  efficient  estimator  and  superior  to  some  other 
approaches  [1,3].  In  [1],  a  single-tone  complex  sinusoid  in 
the  discrete  form  was  treated.  The  frequency  of  the  sinusoid 
was  identified  as  the  one  which  maximizes  the  magnitude  of  the 
Fourier  transform  of  the  discrete  signal.  It  was  obtained  by 
first  searching  coarsely  over  the  DFT  spectrum  and  then 
refining  the  estimate  with  the  cosine  iteration  algorithm. 

This  method  was  extended  to  handle  multiple  tones  in  [2]  where 
the  coarse  search  was  performed  in  a  frequency  space  with 
dimension  equal  to  the  number  of  complex  sinusoids. 

The  ML  method  for  the  complex  data  can  be  used  in  two 
ways  to  treat  real  data  without  any  modification.  One  way  is 
to  regard  each  real  sinusoid  as  a  pair  of  complex  ones  with 
frequencies  equal  in  magnitude  and  opposite  in  sign.  In  doing 
so,  the  dimension  of  the  search  space  is  thus  doubled  and,  as  a 
consequence,  the  computation  cost  would  increase 
significantly.  Furthermore,  the  estimation  precision  would 
also  degrade  because  the  pair  of  frequency  estimates  would  not 
be  necessary  equal  to  each  other  in  magnitude,  especially  in 
the  presence  of  noise  or  in  the  case  of  very  low  frequency 
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The  other  way  is  to  first  generate  an  "analytic" 
(complex-valued)  signal  for  the  given  real -valued  signal  using 
the  Hilbert  transform  [4].  The  difficulty  with  this  approach 
is  that,  for  a  short  discrete  signal,  the  Hilbert  transform 
could  not  produce  a  well-approximated  complex  sinusoid  from 
each  real  one  and  it  would  distort  the  noise  structure. 

In  this  paper  we  present  a  refined  ML  method  for 
frequency  estimation  of  real  sinusoids  in  which  the  dimension 
of  search  space  is  kept  as  the  number  of  real  sinusoids  and  the 
fine  search  is  carried  out  using  a  single-step  interpolation 
which  requires  little  extra  computation.  Simulation  results 
show  that  this  refinement  can  improve  the  estimation  precision 
and  computation  efficiency. 


II. 


THE  METHOD 


A.  Derivation 

We  assume  that  the  observed  signal  {yn,  n=0,..N-1} 
contains  M  real  sinusoids  corrupted  by  additive  noise  wn. 

M 

yn  3  l  amsin  ^<tmn+0m^  +  wn  M> 

m=1 

where  the  real  positive  amplitudes  {amf,  the  phases  {em} 
and  the  frequencies  {u,m}  are  fixed  but  unknown  parameters  and 
{wn}  are  identical  and  independent  Gaussian  distributions 
with  zero  mean  and  variance  aw2 .  Since 

sin  *  -  (e^*-e“^* )/2j 

(1)  can  be  rewritten  in  the  vector-matrix  form  as 

a 

X  *  (P  P')  +  w  (2) 

o' 

where  (yo»* • • »Yn-1 

W  =  (WQ,. • • ,wN_i )t. 

a  ■  (ot , • • om“(am/2j)©^flm 

PNxM  3  1  Pnm  1  •  Pnm*e^n“‘m* 

and  the  superscripts  "t"  and  denote  matrix  transpose  and 
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complex  conjugate  respectively.  By  leting  j?  *(^, )  and  S=(P]P') 
(2)  can  be  further  reduced  to 

y  *  S  £  +  w  (3) 

Note  that  in  the  above  formulation,  each  real  sinusoid  is 
decomposed  into  a  pair  of  complex  ones  whose  frequencies  are 
equal  in  magnitude. 

The  ML  method  estimates  the  unknown  parameters  by 
maximizing  the  likelihood  function  which  is  a  strictly 
increasing  function  of  the  quadratic  form 

Q^tt.,0)  -  -(^-S£)t  (*-S£)  (4) 

whc  rs  (l  *  ( &  ^  f  •  •  f  ^  • 

By  solving  for  £  ,(£  *  (S*S)-1  S*^) ,  and  substituting 

it  back  into  (4)  followed  by  dropping  the  constant  term  “X^X  we 
have 


Q2<£)  =  (xfcS)  (S*S)“1(S*Z) 


(5) 


or 


Q2(±)  =  (X* p  2tp,){ - 


' p *  p  !  p*  p,>r1 


pfc  p  !  pfc  p’> 


P*  x 


Y, 


(6) 


where  the  asterisk  ("*")  denotes  conjugate  transpose. 
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Now  letting 


^Mxl*  N  (p*^) 


AMxN=  N  (P*P) 


BMxM~  N  (p*p,) 


and  inverting  the  partioned  matrix,  we  have 
Q2(w)  =  2N  {d*  R  d  +  Re  [dfcT  d]} 


(7) 


where  R  -  (A-B(A^)“^B' )“^ 
and  T  =  -(At)“1  B'R 

It  is  worthwhile  to  point  out  that  the  component  of  dm  is 
equal  to 


d(w) 


N 


N-1 

I  y 

n=0 


n 


e-jwn 


(8) 


evaluated  at  ustu  and  the  entries  of  A  and  B,  a  and  b  ,  are 
m  mn  mn 

are  equal  to 


P  (u>) 


2 

N 


N-1 

l  e~lmn 
n=0 


sin(Ntt)/2)  -j  (^)w 
N  sin  ( (o/2 )  e  z 


(9) 


evalued  at  ai«w  -<u  and  <u=w  +<u  respectively.  So  Q_(u>)  can  be 
inn  mn  z  _ 
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computed  using  d(u>)  and  p(<o)  only  and  the  frequency  estimate  w 
is  the  w  which  maximizes  this  quantity.  Also  note  that  d(<o)  is 
the  DPT  value  of  the  observed  signal  {yn}  at  u  and  p(w)  is 
the  Dirichlet  kernel  associated  with  the  rectangular  window. 

In  fact  p(w)  also  measures  the  interference  between  two 
spectral  lines  separated  by  u.  A  plot  of  p(cu)  for  N=16r  32  and 
64  is  shown  in  Fig.  1. 

Now  let  us  consider  some  special  cases.  For  the  case 
of  M=1  (  a  single  real  sinusoid) ,  Q2(tt)  can  ^  written  ex¬ 
plicitly  as 


Q2  { w1 )  =  2N{|(dUt)|2-Re[d2(u)l)p'(2<i,1)j}/(1-|p(2uil)|2)  (10) 

Note  that  the  interference  term  p(2mi)  between  the  positive 
and  negative  spectral  lines  of  the  real  sinusoid  appears  in  the 
expression.  In  the  limiting  case  u>i=0,  P=1 
and  Q2=|d(<i)=0)  |  2. 

Suppose  <  h>2  <  wM.  If  is  large  enough  so  that 
the  negative  spectral  lines  associated  with  the  real  sinusoids 
do  not  interfere  with  the  positive  ones  then  B*0  and  Q2  re¬ 
duces  to 


Q2(w)  »  2Nd*A~1d  (11) 

This  is  exactly  the  same  formula  as  Eq.  (51)  of  [2]  for  M 
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complex  sinusoids  except  the  constant  factor  2.  Therefore, 
under  this  condition,  the  frequency  estimates  of  M  real 
sinusoids  can  be  obtained  by  regarding  them  as  M  complex 
sinusoids  with  positive  frequencies. 

From  (11)  it  is  also  clear  that  the  corresponding 
formula  for  one  and  two  complex  sinusoids  are,  after  dropping 
the  factor  2, 


Q2(o)1  )  =  N|d(w1 )  |2 


( 12a) 


and 

Q£(wrw2)  =  N/(  1- 1  p  (fa>2-Ul )  |2  )  x 

{  |d(0i1  )  |  2+|d(  oo2)  |  2-2  Re  ( d  (  a;  1  )  d  1  (  uj2  )  p  (  uj  2~ go  1 ) )  } 


(12b) 


respectively.  Note  that  if  u>2  =  -u^  than  (12b)  reduces  to  (10) 
as  expected. 

B.  Algorithm 

* 

To  find  the  frequency  estimate  w,  using  either  (10) 
or  (11)  we  first  perform,  as  proposed  in  [1-2],  a  finite  search 
in  M  nested  do-loops  for  the  maximum  of  Q2(w.)  over  a 
lattice  when  Q  *  { 2nk/K,k«0 , . . .K-1  and  K>N| .  Empirically 
speaking,  K=2N  or  4N  is  enough  as  long  as  1/K  is  comparable  to 
the  lowest  frequency  in  the  signal  .  The  required  (d ( w )  ,oueQ ) } 
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and  {p(o.),  u.eQ }  during  the  finite  search  can  be  fetched  from 
the  precomputed  (using  the  FFT)  and  prestored  tables.  This 
course  estimate  u.  is  then  refined  by  using  a  single-step 
interpolation  (i.e.  one-iteration  Newton's  method)  to  give  the 
final  estimate 

a  ~  i 

u^u.-H-  q 

where  £  (gradient)  ■  302/3^  and  H( Hessian)  =  are 

evaluted  at  u,  numerically  (based  on  the  central  differences 
(5])  using  values  of  02(0  already  computed  in  the  finite 

A  A 

search.  If  Q2 (fc) <Q2 (u.)  then  ^  is  taken  as  a,.  It  is  found 
that  this  simple  interpolation  can  eliminate  the  lengthy 
iterations  and  the  associated  convergence  problems  and  yet 
provide  a  satisfactory  precision. 

When  M  is  large  it  will  be  difficult  to  carry  out  the 
search  at  one  time  because  of  the  large  search  space  gM. 
Fortunately,  in  many  cases,  M'  (<M)  independent  (non-interfer¬ 
ing)  and  significant  (detectable)  intervals  with  length  o'  <  n 
can  be  identified  from  the  DFT  spectrum  of  the  data, 

{  |d(u. )  |  ,u,eo} ,  each  of  which  contains  one  single  or  a  small 
number  of  blurred  spectral  lines  and  these  intervals  can  be 
treated  individually  in  much  smaller  space  using  either  (7)  or 
(11)  depending  upon  their  distances  from  a.=0. 
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III. 


SIMULATION  RESULTS 


Three  examples  are  used  here  to  test  the  estimation 
precision  and  computation  efficiency  of  the  proposed  ML 
method.  The  first  two  have  been  utilized  in  other  frequency 
estimators  [6,7].  The  computer  simulation  results  for  these 
two  examples  are  compared  with  the  published  data.  In  the  fol¬ 
lowing  the  simulation  conditions  for  each  example  are  specified 
in  terms  of  the  observation  model,  the  number  of  given  samples, 
(N),  the  number  of  DFT  values  used  in  the  search  (K),  the 
number  of  Monte-Carlo  simulations  (MC)  and  the  signal-to-noise 
ratio  which  is  defined  as  SNR  =10  log  (1/2 o£) 

Example  1:  yn  =  sin(.5nn  +.5ir)  +  wn,  N=32 

This  example  was  used  in  [6]  where  the  covariance  and  modified 
covariance  methods  of  the  "maximum  entropy"  (ME)  spectral 
estimator  were  studied.  The  simulation  conditions  are  MC=1000, 
SNR=37  dB  and  K=64.  The  best  performances,  in  terms  of 
estimation  standard  derviations,  of  the  two  ME  methods  as 
extracted  from  Pig.  2  of  [6]  are  .53x10“*  Hz  and  .47x10"*  Hz. 
They  are  outperformed  by  the  ML  method  where  the  standard 
deviation  is  .36x10“4  Hz. 

Example  2:  yn=  sin(.lnn)  +  wn,  N=25 


This  example  was  used  in  [7]  where  the  instrumental  variable 
method  (IVM)  was  proposed  for  estimating  frquencies  of  real 
sinusoids.  The  simulation  conditions  are  MC*200,  K*64  and  SNR 
=22.2,  12.2  and  7  dB.  The  comparison  is  given  in  Table  I. 

Both  Eq.  (10)  and  Eq.  (12b)  are  used.  In  the  latter  case,  the 
real  sinusoid  is  regarded  as  a  pair  of  complex  sinusoids  and 
the  actual  frequency  estimate  is  taken  as  the  average  of  the 
magnitudes  of  the  two  estimated  frequencies  based  on  ( 12b) . 

From  Table  I,  it  is  clear  that  the  ML  method  is 
superior  to  the  IVM  and  its  estimation  precision  approaches  to 
the  Cramer-Rao  bound  ( CRB ) * ( 1 ] .  Although  the  estimation 
precisions  obtained  by  using  (10)  and  (12b)  are  comparable,  the 
required  computing  times  are  quite  different.  On  a  Amdhal  470 
machine,  each  run  takes  . 4050x1 0~2  seconds  using  (10)  and  .1026 
seconds  (25  times  longer)  using  (12b). 

Example  3:  yn  =  sin(.05nn  +  n/4)  +  w 
N=  1 6 

SNR=7 , 5  and  0  dB. 

This  example  is  intended  to  compare  the  performances  of  using 
Eq  (10)  and  Eq.  (12b)  in  the  cases  of  low  SNR's  and  strong 
interference  between  the  positive  and  the  negative  spectral 
lines  of  the  real  sinusoid.  The  interference  is  23%  in  this 

♦The  CRB  is  a  lower  bound  on  the  estimation  variance  of  any 
unbiased  estimator. 

11 


V 


example.  Table  II  gives  the  comparison.  Undoubtedly,  the 
proposed  ML  method  (Eq.  (10))  for  frequency  estimation  of  real 
sinusoid  is  much  more  precise  than  the  one  (Eq.  (12b))  which 
is  adapted  from  the  ML  method  originally  designed  for  complex 
data.  From  Table  II  also  note  that  the  estimation  becomes 
biased  at  very  low  SNR,  e.g.,  0  dB. 


TABLE  I 


Comparison  of  IVM,  MLM  with  Eq.  (10)  and  NLN 
with  (Eq.  12b)  for  frequency  estimation  precision  of  a 
real  sinusoid  in  noise.  The  true  frequency  is  .05  Hz 
Each  entry  in  column  2,3,  and  4  means  meand  ±SD  in 
units  of  Hz.  See  example  2. 


Method 

SNR  dB 

IVM  [6] 

MLM 

Eq.  (10) 

MLM 

Eq.  (12b) 

1/2 

(CRB) 

22.2 

.5013  x  10-1 

±1.897  x  10-3 

.4946  X  10-1 

±.3827  x  10“3 

.4958  x  10-1 

±.3885  x  10-3 

.3244  x  10-3 

12.2 

.5037  x  10_1 

±1.0904  x  10-2 

.4956  x  10~l 

±.1073  x  10-2 

.4972  x  10~l 

±.1089  x  10-2 

.1026  x  10-2 

7.0 

.4819  x  10-1 

±3.0412  x  10“2 

.4962  x  10-1 

±.1941  x  10“2 

.4974  x  10"1 

±.1966  x  10“2 

.1867  x  10“2 

TABLE  II 

Comparison  of  MLM  with  Eq.  (10)  and  MLM  with  Eq.  (12b)  for 
frequency  estimation  precision  of  a  real  sinusoid. 

The  true  frequency  is  .025  Hz.  See  example  3. 


SNR  dB 

method 

7 

5 

0 

Eq.  (10) 

.2313  x  10-1 

±.5199  x  10"2 

.2215  x  10-1 

±.6044  x  10-2 

.3251  x  10"1 

±.5758  x  10*1 

Eq.  (12b) 

.2599  x  10“l 

±.1921  x  10_l 

.3573  x  10_1 

±.4230  x  10"1 

.8304  x  10_l 

±.8218  x  10_1 
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IV. 


CONCLUSIONS 


The  ML  method  [1-2]  for  frequency  estimation  of 
complex  sinusoids  is  extended  to  treat  real  data.  This 
extension  provides  better  estimation  precision  and  computation 
efficiency  than  the  conventional  approaches.  The  detection 
problem(determining  the  number  of  sinusoids  present  in  the 
data)  which  is  not  addressed  here  can  be  handled  easily  based 
on  the  M-ary  generalized  likelihood  ratio  test  [8]  using 
by-products  of  the  ML  estimator. 


| 
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